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pg . A semi-grand-canonical Monte Carlo algorithm is employed in conjunction with the 

bond fluctuation model to investigate the critical properties of an asymmetric binary 
(AB) polymer mixture. By applying the equal peak-weight criterion to the concentration 
distribution, the coexistence curve separating the A-rich and B-rich phases is identified 
as a function of temperature and chemical potential. To locate the critical point of the 
model, the cumulant intersection method is used. The accuracy of this approach for 
determining the critical parameters of fluids is assessed. Attention is then focused on 
the joint distribution function of the critical concentration and energy, which is analysed 

'5^ | using a mixed-field finite-size-scaling theory that takes due account of the lack of symme- 

try between the coexisting phases. The essential Ising character of the binary polymer 
critical point is confirmed by mapping the critical scaling operator distributions onto in- 
dependently known forms appropriate to the 3D Ising universality class. In the process, 
estimates are obtained for the field mixing parameters of the model which are compared 
both with those yielded by a previous method, and with the predictions of a mean field 
calculation. 
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1 Introduction 

The critical point of binary liquid and binary polymer mixtures, has been a subject of abiding 
interest to experimentalists and theorists alike for many years now. It is now well established 
that the critical point properties of binary liquid mixtures fall into the Ising universality class 
(the default for systems with short ranged interactions and a scalar order parameter) [p]]. 
Recent experimental studies also suggest that the same is true for polymer mixtures f^j-[|TT|1. 
However, since Ising-like critical behaviour is only apparent when the correlation length far 
exceeds the polymer radius of gyration £ 3> (Rg), the Ising regime in polymer mixtures 
is confined (for all but the shortest chain lengths) to a very narrow temperature range near 
the critical point. Outside this range a crossover to mean-field type behaviour is seen. The 
extent of the Ising region is predicted to narrow with increasing molecular weight in a manner 
governed by the Ginsburg criterion [12|, disappearing entirely in the limit of infinite molecular 



weight. Although experimental studies of mixtures with differing molecular weights appear to 
confirm qualitatively this behaviour ||, there are severe problems in understanding the scaling 
of the so-called "Ginsburg number" (which marks the centre of the crossover region and is 
empirically extracted from the experimental data |l3j) with molecular weight 

Computer simulation potentially offers an additional source of physical insight into polymer 
critical behaviour, complementing that available from theory and experiment. Unfortunately, 
simulations of binary polymer mixtures are considerable more exacting in computational terms 
than those of simple liquid or magnetic systems. The difficulties stem from the problems of deal- 
ing with the extended physical structure of polymers. In conventional canonical simulations, 
this gives rise to extremely slow polymer diffusion rates, manifest in protracted correlation 
times [14], [jl|. Moreover, the canonical ensemble does not easily permit a satisfactory treat- 



ment of concentration fluctuations, which are an essential feature of the near-critical region in 
polymer mixture. In this latter regard, semi-grand-canonical ensemble (SGCE) Monte Carlo 
schemes are potentially more attractive than their canonical counterparts. In SGCE schemes 
one attempts to exchange a polymer of species A for one of species B or vice-versa, thereby 
permitting the concentration of the two species to fluctuate. Owing, however, to excluded 
volume restrictions, the acceptance rate for such exchanges is in general prohibitively small, 
except in the restricted case of symmetric polymer mixtures, where the molecular weights of 
the two coexisting species are identical (Na = Nb)- All previous simulation work has therefore 
focussed on these symmetric systems, mapping the phase diagram as a function of chain length 



and confirming the Ising character of the critical point |16], [T7], [L8| . Tentative evidence for a 

crossover from Ising to mean field behaviour away from the critical point was also obtained 

HJ. Hitherto however, no simulation studies of asymmetric polymer mixtures (Na 7^ Nb) 



have been reported. 

Recently one of us has developed a new type of SGCE Monte Carlo method that amelio- 
rates somewhat the computational difficulties of dealing with asymmetric polymer mixtures 
[POU . The method, which is described briefly in section |3.1| , permits the study of mixtures of 
polymer species of molecular weight Na and Nb = kNA, with k — 2,3, 4- ••. In this paper 
we shall employ the new method to investigate the critical behaviour of such an asymmet- 
ric polymer mixture. In particular we shall focus on those aspects of the critical behaviour of 
asymmetric mixtures that differ from those of symmetric mixtures. These difference are rooted 
in the so called 'field mixing' phenomenon, which manifests the basic lack of energetic (Ising) 
symmetry between the coexisting phases of all realistic fluid systems. Although it is expected 
to have no bearing on the universal properties of fluids, field mixing does engender certain 
non-universal effects in near-critical fluids. The most celebrated of these is a weak energy-like 
critical singularity in the coexistence diameter [F2I], [33], the existence of which constitutes a 



failure for the 'law of rectilinear diameter'. As we shall demonstrate however, field mixing 
has a far more legible signature in the interplay of the near-critical energy and concentration 
fluctuations, which are directly accessible to computer simulation. 

In computer simulation of critical phenomena, finite-size-scaling (FSS) techniques are of 
great utility in allowing one to extract asymptotic data from simulations of finite size p3| . One 



particularly useful tool in this context is the order parameter distribution function [24, 25, 26 



Simulation studies of magnetic systems such as the Ising [25] and 4 models [27|, demonstrate 
that the critical point form of the order parameter distribution function constitutes a useful 
hallmark of a university class. Recently however, FSS techniques have been extended to fluids 
by incorporating field mixing effects p8| , |29j. The resulting mixed-field FSS theory has been 
successfully deployed in Monte Carlo studies of critical phenomena in the 2D Lennard- Jones 
fluid |5?| and the 2D asymmetric lattice gas model [[!tj . 



The present work extends this programme of field mixing studies to 3D complex fluids with 
an investigation of an asymmetric polymer mixture. The principal features of our study are as 
follows. We begin by studying the order parameter (concentration) distribution as a function 
of temperature and chemical potential. The measured distribution is used in conjunction 
with the equal peak weight criterion to obtain the coexistence curve of the model. Owing to 
the presence of field mixing contributions to the concentration distribution, the equal weight 
criterion is found to break down near the fluid critical point. Its use to locate the coexistence 
curve and critical concentration therefore results in errors, the magnitude of which we gauge 
using scaling arguments. The field mixing component of the critical concentration distribution 
is then isolated and used to obtain estimates for the field mixing parameters of the model. 
These estimates are compared with the results of a mean field calculation. 

We then turn our attention to the finite-size-scaling behaviour of the critical scaling op- 
erator distributions. This approach generalises that of previous field mixing studies which 
concentrated largely on the field mixing contribution to the order parameter distribution func- 
tion. We show that for certain choices of the non-universal critical parameters — the critical 
temperature, chemical potential and the two field mixing parameters — these operator distri- 
butions can be mapped into close correspondence with independently known universal forms 
representative of the Ising universality class. This data collapse serves two purposes. Firstly, it 
acts as a powerful means for accurately determining the critical point and field mixing param- 
eters of model fluid systems. Secondly and more generally, it serves to clarify the sense of the 
universality linking the critical polymer mixture with the critical Ising magnet. We compare 
the ease and accuracy with which the critical parameters can be determined from the data 
collapse of the operator distributions, with that possible from studies of the order parameter 
distribution alone. It is argued that for critical fluids the study of the scaling operator distri- 
butions represent the natural extension of the order parameter distribution analysis developed 
for models of the Ising symmetry. 

2 Background 

In this section we review and extend the mixed-field finite-size-scaling theory, placing it within 
the context of the present study. 

The system we consider comprises a mixture of two polymer species which we denote A and 
B, having lengths N^ and Nb monomers respectively. The configurational energy $ (which we 
express in units of /c#T) resides in the intra-and inter molecular pairwise interactions between 
monomers of the polymer chains: 



$({r})= £ wdn-rjl) (2.1) 

where M = uaNa + n^Ns with ua and ng the number of A and B type polymers respectively. 
A/" is therefore the total number of monomers (of either species), which in the present study 
is maintained strictly constant. The inter-monomer potential v is assigned a square- well form 

v(r) = -e r<r m (2.2) 

v(r) =0 r > r m 

where e is the well depth and r m denotes the maximum range of the potential. In accordance 



with previous studies of symmetric polymer mixtures [jR], pTTH , we assign e = €aa = ^bb = 
-e AB > 0. 

The independent model parameters at our disposal are the chemical potential difference per 
monomer between the two species A/i = \ia — Hb, and the well depth e (both in units of k B T). 
These quantities serve to control the observables of interest, namely the energy density u and 
the monomer concentrations <pA and </>#. Since the overall monomer density <p^ = 4>a + <Pb is 
fixed, however, it is sufficient to consider only one concentration variable 0, which we take as 
the concentration of A-type monomers: 

= <\> A = L~ d n A N A (2.3) 

The dimensionless energy density is defined as: 

u = L- d e _1 $({r}) (2.4) 

with d — 3 in the simulations to be chronicled below. 

The critical point of the model is located by critical values of the reduced chemical potential 
difference A/i c and reduced well-depth e c . Deviations of e and A/x from their critical values 
control the sizes of the two relevant scaling field that characterise the critical behaviour. In 
the absence of the special symmetry prevailing in the Ising model, one finds that the relevant 
scaling fields comprise (asymptotically) linear combinations of the well-depth and chemical 
potential difference |2~T| ]: 

r = e c — e + s(A/i — A/x c ) h = A/i — A/i c + r(e c — e) (2-5) 

where r is the thermal scaling field and h is the ordering scaling field. The parameters s 
and r are system-specific quantities controlling the degree of field mixing. In particular r is 
identifiable as the limiting critical gradient of the coexistence curve in the space of A/i and e. 
The role of s is somewhat less tangible; it controls the degree to which the chemical potential 
features in the thermal scaling field, manifest in the widely observed critical singularity of the 
coexistence curve diameter of fluids |1|, [22|, jjflfl . 



Conjugate to the two relevant scaling fields are scaling operators £ and J\A, which comprise 
linear combinations of the concentration and energy density [ER 123] : 



M = j^gf [(p - su] 6 = j^gr \ u ~ r( t ) \ ( 2 -6) 

The operator A4 (which is conjugate to the ordering field h) is termed the ordering operator, 
while S (conjugate to the thermal field) is termed the energy-like operator. In the special case 



of models of the Ising symmetry, (for which s = r = 0), M. is simply the magnetisation while 
£ is the energy density. 

Near criticality, and in the limit of large system size L, the probability distributions Pl{M) 
and Pl(£) of the operators M. and £ are expected to be describable by finite-size-scaling 
relations having the form [^, ^| : 

p L (M) ~ a M - l L d ^ M p M {a M - x L d - XM bM,a M L XM h,a e L X£ T) (2.7a) 

p L (£) ~ a £ - 1 L d - X£ p £ (a £ - l L d - X£ 5£,a M L XM h,a £ L X£ T) (2.7b) 

where 5.M = M. — M. c and 5£ = £ — £ c . The functions p^ and p £ are predicted to be 
universal, modulo the choice of boundary conditions and the system-specific scale-factors cl_m 
and as of the two relevant fields, whose scaling indices are \m = d — /3/u and A^ = \jv 
respectively. Precisely at criticality (h = r — 0) equations p.7a| and |2.7b| imply 



p L {M) ~ a M - x L^ v f M {a M - x L^ v bM) (2.8a) 

p L (£) ~ o e - 1 L (1 - a)/,, ^(afi- 1 L (1 - a)/, '55) (2.8b) 



where 



p^(x) =p^(x,y = 0,z = 0) #Ka;)=|fg(a;,3/ = 0,z = 0) (2.9) 

are functions describing the universal and statistically scale-invariant fluctuation spectra of 
the scaling operators, characteristic of the critical point. 

The claim that the binary polymer critical point belongs to the Ising universality class is 
expressed in its fullest form by the requirement that the critical distribution of the fluid scaling 
operators Pl(M) and Pl{£) match quantitatively their respective counterparts — the magneti- 
sation and energy distributions — in the canonical ensemble of the critical Ising magnet. As we 
shall demonstrate, these mappings also permit a straightforward and accurate determination 
of the values of the field mixing parameters s and r of the model. 

An alternative route to obtaining estimates of the field mixing parameters is via the field 
mixing correction to the order parameter (i.e. concentration) distribution pl{4>)- At criticality, 
this distribution takes the form [E9l BUI: 



p L (<t>) ~ a M ~ X ^ lh 
where 



p* M {x) - sa £ a M - l L-^- a -^^- (pU(x)u>*(x)) + 0(s 2 [ 



\x=a M - x L^l"\4>-4> c \ 
(2.10) 



u*{x) = a^L'-^iiu^)) -u c - r(0 - C )] + O(s) (2.11) 



which we term the energy function, is a universal function characterising the critical energy- 
like operator. Equation [2.10| states that to leading order in the field mixing parameter s, 



the order parameter distribution is a sum of two distinct universal components. The first of 



these, p* M {x) } is simply the universal ordering operator distribution featuring in equation |2.8a . 



The second, J| (p% i (x)uj*(x)), is a function characterising the mixing of the critical energy- 
like operator into the order parameter distribution. This field mixing term is down on the 
first term by a factor j_ i -\ x -° l ~p)I v anc [ therefore represents a correction to the large L lim- 
iting behaviour. Given further the symmetries of uj(x) and p* M (x), both of which are even 
(symmetric) in the scaling variable x | |29|j , the field mixing correction is the leading antisym- 



metric contribution to the concentration distribution. Accordingly, it can be isolated from 



measurements of the critical concentration distribution simply by antisymmetrising around 
4>c = (0)c- The values of s and r are then obtainable by matching the measured critical func- 
tion — s-^T {pl(4>) KM0)) ~~ u c — r (0 — 0c)]} to the measured antisymmetric component of the 



critical concentration distribution [30]. In the simulations described below we shall compare 
the values of s and r obtained by this method with those obtained by matching the fluid 
operator distributions to their Ising equivalent forms. 

Finally in this section, we turn to a consideration of the finite-size-scaling behaviour of the 
energy distribution function. In contrast to the situation for the order parameter distribution 
described above, it transpires that field mixing radically alters the limiting form of the critical 
energy density distribution. To substantiate this claim we reexpress u in terms of the scaling 



operators. Appealing to equation [2.6| , one finds 



u = S + rM (2.12) 

so that the critical energy density distribution is 

Pl(u) = Pl(S + rM) (2.13) 



Now the structure of the scaling forms |2.8aj and |2.8b| show that the typical size of the fluc- 



tuations in the energy-like operator will vary with system size like 58 ~ L-( x - a )l v ^ while the 
typical size of the fluctuations in the ordering operator vary like 5Ai ~ L~^l v . Given that 
in general a < (3, it follows that asymptotically, the contribution of £ to the argument on the 
right hand side of equation |2.13| can be neglected, so that 



p L (u) ~ PL (rM) ~ a M - l rLV»p* M (a M - l rLPI v 5M) (2.14) 

We conclude then that for sufficiently large L, the distribution of the fluid critical energy 
density has the same functional form as the distribution of the critical ordering operator 
Pm . Given further that p* M possesses a symmetric double-peaked form, while the critical 
energy distribution in the Ising model Pl{u) = p£ possesses an asymmetric single-peaked form, 
the profound influence of field mixing on the critical energy distribution of fluids should be 
apparent. 

3 Monte Carlo studies 

3.1 Algorithmic and computational aspects 

The polymer model studied in this paper is the bond- fluctuation model (BFM). The BFM is a 
coarse-grained lattice-based model that combines computational tractability with the impor- 
tant qualitative features of real polymers: monomer excluded volume, monomer connectivity 
and short range interactions. Within the framework of the model, each monomer occupies a 
whole unit cell of a 3D periodic simple cubic lattice. Neighbouring monomers along the poly- 
mer chains are connected via one of 108 possible bond vectors. These bond vectors provide 
for a total of 5 different bond lengths and 87 different bond angles. Thermal interactions are 
catered for by a short range inter-monomer potential. Further details concerning the model 
can be found in reference [32]. 

The system we have studied comprises two polymer species A and B having chain lengths 
Na and Nb, with Nb = kN^- The SGCE scheme whereby polymers of type A are transformed 
into polymers of type B (or vice versa) is described in appendix 1 and in reference p0| |, but in 



general terms operates as follows. Using a Metropolis algorithm, an A-type polymer is formed 



simply by cutting a B-type polymer into k equal segments. Conversely, a B-type polymer is 
manufactured by connecting together the ends of k A-type polymers. This latter operation 
is, of course, subject to condition that the connected ends satisfy the bond restrictions of the 
BFM. Consequently it represents the limiting factor for the efficiency of the method, since for 
large values of k and Na, the probability that k polymer ends simultaneously satisfy the bond 
restrictions becomes prohibitively small. The acceptance rate for SGCE moves is also further 
reduced by factors necessary to ensure that detailed balance is satisfied. In view of this we 
have chosen k = 3, Na = 10 for the simulations described below, resulting in an acceptance 
rate for SGCE moves of approximately 14%. 

In addition to the compositional fluctuations associated with SGCE moves, it is also nec- 
essary to relax the polymer configurations at constant composition. This is facilitated by 
monomer moves which can be either of the local displacement form, or of the reptation ('slith- 
ering snake') variety ||. These moves were employed in conjunction with SGCE moves, in the 
following ratios : 

local displacement : reptation : semi-grandcanonical = 4 : 12 : 1 

the choice of which was found empirically to relax the configurational and compositional modes 
of the system on approximately equal time scales. 

In the course of the simulations, a total of five system sizes were studied having linear 
extent L = 32, 40, 50, 64 and 80. An overall monomer filling fraction of 8<f)^f = 0.5 was chosen, 



representative of a dense polymer melt [H|. Here the factor of 8 constitutes the monomeric 
volume, each monomer occupying 8 lattice sites. The cutoff range of the inter-monomeric 
square well potential was set at r m = ^/6 (in units of the lattice spacing), a choice which 
ensures that the first peak of the correlation function is encompassed within the range of the 
potential. The observables sampled in the simulations were the A-monomer concentration 



and the energy density u (cf. equations |273| and |2T4j ). The joint distribution p L (4>,u) of these 



quantities was accumulated in the form of a histogram, with successive samples being separated 
by 12.5 semi-grand canonical sweeps in order to reduce correlations. The final histograms for 
the lattice sizes L = 40 and 64 each comprised some 5 x 10 5 entries. Assistance in exploring 
the phase space of the model, was provided by means of the multi-histogram reweighting 
technique |33], [34]. This technique allows one to generate estimated histograms pl(4>,u) for 



values of the control parameters e and A/x other than those at which the simulations were 
actually performed. Such extrapolations are generally very reliable in the neighbourhood of 
the critical point, due to the large critical fluctuations [[33)1 . 

3.2 The coexistence curve and critical limit 

In general for fluid systems, the coexistence curve is not known a-priori and must therefore be 
identified empirically as a prelude to locating the critical point itself. One computational crite- 
rion that can be used to effect this identification is the so-called 'equal weight criterion' for the 
order parameter (concentration) distribution function Pl{4>) = J dupL(<p,u) ||35| . Precisely on 
coexistence and for temperatures well below criticality, Pl{<P) will comprise two well-separated 
gaussian peaks of equal weight, but unequal heights and widths. The centres of these peaks 
identify the concentrations of the coexisting A-rich and A-poor phases. For a given subcritical 
well-depth e, the coexistence value for the chemical potential difference Afi cx can therefore be 
obtained by adjusting A/i until the concentration distribution satisfies the condition: 

r<t>M 
p L ((j), e,Afx cx ) d(f)= / Pl(0,£, A/i ca: ) dcf) (3.1) 



where 0* is a parameter defining the boundary between the two peaks. 

Well below criticality, the value of A/x cx obtained from the equal weight criterion is in- 
sensitive to the choice of 0*, provided it is taken to lie approximately midway between the 
peaks and well away from the tails. As criticality is approached however, the tails of the 
two peaks progressively overlap making it impossible to unambiguously define a peak in the 
manner expressed by equation |3J]. For models of the Ising symmetry, for which the peaks are 
symmetric about the coexistence concentration <p cx , the correct value of Afi cx can nevertheless 
be obtained by choosing 0* = (0) in equation JO] , In near-critical fluids, however, the imposed 
equal weight rule forces a shift in the chemical potential away from its coexistence value in or- 
der to compensate for the presence of the field mixing component. Only in the limit as L — > oo 
(where the field mixing component dies away), will the critical order parameter distribution 
be symmetric allowing one to choose 0* = ((f)) and still obtain the correct coexistence chemical 
potential. Thus for finite-size systems, use of the equal weight criterion is expected to lead 
to errors in the determination of A/z cx near the critical point. Although this error is much 
smaller than the uncertainty in the location of the critical point along the coexistence curve 
(see below), it can lead to significant errors in estimates of the critical concentration C . 

To quantify the error in C it is necessary to match the magnitude of the field mixing 
component of the concentration distribution w(Spl), to the magnitude of the peak weight 
asymmetry w'(8n) associated with small departures 8/j, = A/x — Ap cx from coexistence: 

w(8p L ) = v/(6fi) (3.2) 



Now from equation |2.10| 



w(8p L ) « f M d<f> 5p L {<P) ~ l^ 1 -"-®/" (3.3) 



while 

w'(8p) « [*" d<j> %^ tf/i oc L«^>/" 8p (3.4) 

J<j>* oAfj, 

It follows that the error in A/i c:r varies with system size like: 

6fi oc L- {1 - a+ ^ )/u (3.5) 

Accordingly the error in the critical concentration (obtained as the first moment of the con- 
centration distribution) varies with system size like 

<50 c = x(L)6fM = U lv b\i ~ L- {1 - a)lu (3.6) 

Note also that an analogous treatment of the shift in the auxiliary variable <fi* leads to the 
same L-dependence. 

Measurements of the concentration distribution were performed in conjunction with the 
equal weight criterion, to locate the coexistence curve as a function of well depth and chemical 
potential. The results are shown is figure |l[ Since in finite-size systems, the order parameter 
distribution exhibits a double peaked structure even above the critical temperature, the data 
shown also represent the analytic continuation of the true coexistence curve that persists in 
finite-size systems [p9[ . To determine the position of the critical point on this line of pseudo 



coexistence, the cumulant intersection method was employed. The fourth order cumulant ratio 
Gl is a dimensionless quantity that characterises the form of a function. It is defined as 

°l = 1 " ^f% (3-7) 



where m? and m 4 are the second and fourth moments respectively of the order parameter 
m = — (0). To the extent that field mixing corrections can be neglected, the critical 
order parameter distribution function is expected to assume a universal scale invariant form. 
Accordingly, when plotted as a function of e, the coexistence values of Gl for different system 



sizes are expected to intersect at the critical well depth e c [p6| . This method is particularly 
attractive for locating the critical point in fluid systems because the even moments of the order 
parameter distribution are insensitive to the antisymmetric (odd) field mixing contribution. 
Figure |2] displays the results of performing this cumulant analysis. A well-defined intersection 
point occurs for a value Gl = 0.47, in accord with previously published values for the 3D Ising 
universality class ||36|| . The corresponding estimates for the critical well depth and critical 
chemical potential are 

e c = 0.02756(15) A/x c = 0.003603(15) 

It is important in this context, that a distinction be drawn between the errors on the location 
of the critical point, and the error with which the coexistence curve can be determined. The 
uncertainty in the position of the critical point along the coexistence curve, as determined from 
the cumulant intersection method, is in general considerably greater than the uncertainty in 
the location of the coexistence curve itself. This is because the order parameter distribution 
function is much more sensitive to small deviations off coexistence (due to finite e — e cx or finite 
Afj, — A/x cx ) than it is for deviations along the coexistence curve, (e and A/x tuned together 
to maintain equal weights). In the present case, we find that the errors on A/x c and e c are 
approximately 10 times those of the coexistence values e cx and Afi cx near the critical point. 

The concentration distribution function at the assigned value of e c and the corresponding 
value of AfjLcx, (determined according to the equal weight rule with <fi* =< <p >), is shown 
in figure [3] for the L = 40 and L = 64 system sizes. Also shown in the figure is the critical 
magnetisation distribution function of the 3D Ising model obtained in a separate study [pTj] . 



Clearly the L = 40 and L = 64 data differ from one another and from the limiting Ising 
form. These discrepancies manifest both the pure antisymmetric field mixing component of 
the true (finite-size) critical concentration distribution, and small departures from coexistence 
associated with the inability of the equal weight rule to correctly identify the coexistence 
chemical potential. To extract the infinite-volume value of <p c from the finite-size data, it is 
therefore necessary to extrapolate to the thermodynamic limit. To this end, and in accordance 



with equation 3.6 , we have plotted CX (L), representing the first moment of the concentration 
distribution determined according to equal weight criterion at the assigned value of e c , against 
jJx-ol)/v _ This extrapolation (figure |j) yields the infinite-volume estimate: 

C = 0.03813(19) 

corresponding to a reduced A-monomer density (j) c /(j)j^ = 0.610(3). The finite-size shift in the 
value of 4>cx(L) is of order 2%. 

We turn next to the determination of the field mixing parameters r and s. The value of r 
represents the limiting critical gradient of the coexistence curve which, to a good approxima- 
tion, can be simply read off from figure p] with the result r = —0.97(3). Alternatively (and as 
detailed in ||29||) r may be obtained as the gradient of the line tangent to the measured critical 



energy function (equation |2.11| ) at <p = <p c . Carrying out this procedure yields r = —1.04(6). 

The procedure for extracting the value of the field mixing parameter s from the concentra- 
tion distribution is rather more involved, and has been described in detail elsewhere [^, [3(] 
The basic strategy is to choose s such as to satisfy 

d 
5Pl{4>) = s— (p L (0) [(m(0)) -u c - r{(j> - C )]} 



where 8pl{4>) is the measured antisymmetric field mixing component of the critical concen- 
tration distribution |30|, obtained by antisymmetrising the concentration distribution about 



4> C (L) and subtracting additional corrections associated with small departures from coexistence 
resulting from the failure of the equal weight rule. Carrying out this procedure for the L = 40 
and L = 64 critical concentration distributions yields the field mixing components shown in 
figure |5[ The associated estimate for s is 0.06(1). Also shown in figure |5] (solid line) is the 
predicted universal form of the 3D order parameter field mixing correction —£- (pj^(x)u)*(x)) 



(cf. equation |2.10| ) obtained from independent Ising model studies [[J7[|. Clearly the measured 
functional form of the field mixing correction is in reasonable accord with the universal pre- 
diction. We attribute the residual discrepancies to field mixing contributions of order s 2 or 



higher, not included in equation |2.10| . The directions of the two relevant scaling fields h and r 



corresponding to the measured values of s and r are indicated on figure [TJ. 

3.3 The critical limit revisited : Scaling operator distributions 

In this subsection we consider an alternative method for locating the critical point and deter- 
mining the field mixing parameters s and r , that circumvents some of the difficulties associated 
with using the order parameter distribution alone. 

The method focuses on the ordering and energy-like scaling operator distributions Pl{M) 
and Pl{£) (cf- equations |2.8a| , 2.8b ), which are obtained from the joint distribution of the 



energy density and concentration pi((j),u) as 

M M) =P J^) M£)=PL (lZll) (3.8) 

VI — sr J \\ — sr J 

Precisely at criticality, Pl{-M) is expected to match the universal fixed point function p* M . 
This suggests that the critical parameters can be readily located by simultaneously adjusting 
s,e and A/x until pl{[4> — su]/{l — sr)) matches p* M (modulo the choice of non-universal scale 
parameters implicit in the definition of the scaling variable). The results of performing this 
procedure for pl(-M) are displayed in figure || for the L = 40 and L = 64 system sizes. The 
quality of the data collapse lends substantial support to the contention that the binary polymer 
critical point does indeed belong to the Ising universality class. The corresponding values of 
the critical parameters e, A/x c and s are 

e c = 0.02756(15) Afi c = 0.003603(15) s = 0.06(1) (3.9) 

in good agreement with those determined previously. We note however that the present method 
permits the determination of s without the need to isolate the field mixing component of the 
concentration distribution, a procedure that is somewhat cumbersome and which is anyway 
only accurate to leading order in s |3(J . 

The value of the field mixing parameter r , is intimately associated with the critical energy 
distribution Pl{u), the form of which is shown in figure [7] for both the L = 40 and L = 
64 system sizes. The corresponding mapping of the scaling operator distribution function 
Pl{[u — ?"0]/(l — sr )) onto the universal energy distribution of the 3D Ising model, p\ is shown 
in figure |8]. Again the agreement with the universal form is gratifying, although there are small 
discrepancies which we attribute to corrections to scaling. The data collapse implies a value 
r = —1.00(3), which also agrees to within error with the value obtained previously. 

With regard to the critical energy distributions of figure 0, we note that the distributions 
are not at all reminiscent of the critical Ising energy distribution of figure [5[ Neither, however, 
are they similar to p* M (cf. figure |6|), which it was claimed they match in the limit L — > oo (cf. 



section |2[). This discrepancy implies that the system size is still too small to reveal the asymp- 
totic behaviour Nevertheless the data do afford a test of the approach to the limiting regime, 
via the FSS behaviour of the variance of the energy distribution. Recalling equation |2.14| , 
we anticipate that this variance exhibits the same FSS behaviour as the Ising susceptibility, 
namely: 

L d ((u 2 )-u 2 c )~V /u (3.10) 

By contrast, the variance of the scaling operator S is expected to display the FSS behaviour 
of the Ising specific heat: 

L d ((£ 2 )-£ 2 c )~L a /\ (3.11) 

Figure |^ shows the measured system size dependence of these two quantities at criticality. Also 
shown is the scaled variance of the ordering operator L d ((j\4 2 ) — M. 2 ) ~ L 1 ^ . Straight lines 
of the form L 7//iy and L a ^ u , (indicative of the FSS behaviour of the Ising susceptibility and 
specific heat respectively) have also been superimposed on the data. Clearly for large L, the 
scaling behaviour of the variance of the energy distribution does indeed appear to approach 
that of the ordering operator distribution. 

4 Mean field calculations 

In this section we derive approximate formulae for the values of the field mixing parameters s 
and r on the basis of a mean field calculation. 

Within the well-known Flory-Huggins theory of polymer mixtures, the mean-field equation 
of state takes the form: 

A// = -L ln(p) + -1- ln(l -p)- 2ze(2p -1)+C (4.1) 

In this equation, z ~ 2.7 is the effective monomer coordination number, whose value we have 
obtained from the measured pair correlation function, p = <p/<f)j^ is the density of A-type 
monomers and the constant C is the entropy density difference of the pure phases, which is 
independent of temperature and composition. In what follows we reexpress p by the concen- 
tration 0. 

The critical point is defined by the condition: 

A A „ c = !_ A ^ = o (4.2) 

where A/j, c = A/z(0 c , e c ). This relation can be used to determine the critical concentration and 
critical well-depth, for which one finds 

0C * and I = g/ A NAN ° (4-3) 



(f>N 1 + 1/Vfc e c (VN^+^N B 



,2 



Below the critical point a first order phase transition occurs between the A-rich and A-poor 
phases. To determine the location of the phase boundary we employ a Landau expansion of 
the equation of state in terms of the small parameters 5<f) = <ft — <ft c and 5e = e — e c 



-c- 



Afi = Afi c + r'5e - a5e5<j) + b5<f) s + c<50 4 H (4.4) 



where the expansion coefficients take the form 

, 0c 4^ (l + Vk) A (k-l)(l + Vk) 4 . . 



The phase boundary itself is specified by the binodal condition 

ff+ d0 Au(<f>, e) 
A^ 6 = A/*(0 +l 6) = A/i(0_, 6) = J0 - ^ ^ (4.6) 

(P+-(P- 

where 0_ and + denote the concentration of A monomers in the A-poor phase and A-rich 
phases respectively. Thus to leading order in e, the phase boundary is given by : 

A/i ca .(e) = Afj, c + r'5e c + --- (4.7) 

Consequently we can identify the expansion coefficient r' with the field mixing parameter r 
(c.f. equation |2.5| ) that controls the limiting critical gradient of the coexistence curve in the 
space of Afj, and e. Substituting for A/x c and e c in equation fl7| and setting k = 3, we find 
r = —1.45, in order-of-magnitude agreement with the FSS analysis of the simulation data. 

In order to calculate the value of the field mixing parameter s, it is necessary to obtain 
the concentration and energy densities of the coexisting phases near the critical point. The 
concentration of A-type monomers in each phase is given by 



, a $ e 2ac5e . A . 

5<P ± = ± - C = ±y wT + --- (4-8) 

so that the variation of the order parameter along the coexistence curve is: 

/tj.\ [0+ + 0- , ] 2ac5e 

<*> = {— — *-) = ~sr (49) 

A similar calculation for the energy density yields: 

<«(*)> = "^ (z. + z(2p - I) 2 ) = -^ L + z(2 A _ i)^ + rH _ 2£.^ (410) 



where z s is the coordination number of the intra-chain thermal interactions. The variation of 
the energy density along the coexistence curve then follows as: 

u((f)+) + u{(j)J) 2za ( c(j)x\ x ( . 

(6u) = M (0 C ) = -_ I 1 + r— 1 *e + • • • (4.11) 



Now since (1 — rs)(vM) = (6(f)) — s(Su) vanishes along the coexistence line, equations fO| and 
4.11| yield the following estimate for the field mixing parameter s: 



s = iM = c ^ = 3(fc-l) , A 12 , 

(Su) 5zb(l + r^) 2(W*(l + rfg$) 

Thus within the mean field framework, the field mixing parameter s is controlled by the 
ratio of the fifth and fourth order coefficient of the Landau expansion of the free energy. For 
the present case (k = 3) equation |4.12| yields s = 0.070, in good agreement with the value 



obtained from the FSS analysis of the simulation data. It is also similar in magnitude to the 



values of s measured for the 2D Lennard- Jones fluid 29 and 2D asymmetric lattice gas model 



|30|j . The sign of the product rs differs however from that found at the liquid- vapour critical 
point. In the present context this product is given by 



' >x " x """ (4.13) 



1 - rs 10 VN A N B 

However an analogous treatment of the van der Waals fluid predicts a positive sign rs, in 
agreement with that found at the liquid vapour critical point [^, [30[] . 

5 Concluding remarks 

In summary we have employed a semi-grand-canonical Monte Carlo algorithm to explore the 
critical point behaviour of a binary polymer mixture. The near-critical concentration and scal- 
ing operator distributions have been analysed within the framework of a mixed-field finite-size 
scaling theory. The scaling operator distributions were found to match independently known 
universal forms, thereby confirming the essential Ising character of the binary polymer critical 
point. Interestingly, this universal behaviour sets in on remarkably short length scales, being 
already evident in systems of linear extent L = 32, containing only an average of approximately 
100 polymers. 

Regarding the specific computational issues raised by our study, we find that the concen- 
tration distribution can be employed in conjunction with the cumulant intersection method 
and the equal weight rule to obtain a rather accurate estimate for the critical temperature 
and chemical potential. The accuracy of this estimate is not adversely affected by the anti- 
symmetric (odd) field mixing contribution to the order parameter distribution, since only even 
moments of the distribution feature in the cumulant ratio. Unfortunately, the method can 
lead to significant errors in estimates of the critical concentration C , which are sensitive to 
the magnitude of the field mixing contribution. The infinite- volume value of <fi c must therefore 
be estimated by extrapolating the finite-size data to the thermodynamic limit (where the field 
mixing component vanishes). Estimates of the field mixing parameters s and r can also be 
extracted from the field mixing component of the order parameter distribution, although in 
practice we find that they can be determined more accurately and straightforwardly from the 
data collapse of the scaling operators onto their universal fixed point forms. 

In addition to clarifying the universal aspects of the binary polymer critical point, the 
results of this study also serve more generally to underline the crucial role of field mixing in 
the behaviour of critical fluids. This is exhibited most strikingly in the form of the critical 
energy distribution, which in contrast to models of the Ising symmetry, is doubly peaked with 
variance controlled by the Ising susceptibility exponent. Clearly therefore close attention must 
be paid to field mixing effects if one wishes to perform a comprehensive simulation study of 
critical fluids. In this regard, the scaling operator distributions are likely to prove themselves 
of considerable utility in future simulation studies. These operator distributions represent the 
natural extension to fluids of the order parameter distribution analysis deployed so successfully 
in critical phenomena studies of (Ising) magnetic systems. Provided therefore that one works 
within an ensemble that affords adequate sampling of the near-critical fluctuations, use of the 
operator distribution functions should also permit detailed studies of fluid critical behaviour. 
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Appendix A 

In this appendix we give a brief description of our semi-grandcanonical algorithm for polymer 
mixtures with chain lengths Nb = kN^,k = 3. A more detailed presentation of the method 
can be found in reference | j2U[| . 



As illustrated in figure [10] the SGCE Monte-Carlo moves consists in either joining together 
k A-polymers to form a B chain, or alternatively, cutting a B-polymer into k equal segments, 
each of which is then an A chain. The B chain formation step is attempted with probability 
tiaI^a + u b), and proceeds as follows. First one starts by choosing a random A-polymer 
end (a given end is selected with probability 1/2ua)- One then determines the number v\ of 
neighbouring A-ends that satisfy the bonding constraints. Of the V\ possible ends to which a 
bond might be formed, one is chosen randomly and the ends connected together. In the same 
way one computes the number z/ 2 of possible bonding partners for the remaining end of the 
second A-polymer and makes a connection if possible. Thus the proposition probability for 
B polymer formation is given by: P^Zb' = l/2{ n A + ^b)^i^2- Finally the move is accepted 
with probability 

Pj%_ BI = mm(l,exp(-/3AE(kA -> B') - N B Afj) 

The formation procedure for A chains simply involves cutting a B-chain into 3 equal parts, 
a procedure which is attempted with probability UB'jijiA + u b)- One end of a B chain is 
chosen with probability 1/2tir, and the cutting procedure starts from this end, leading to a 
proposition probability: P^'ZkA = l/2{ n A + wg). For the acceptance, one determines the 
possible number of bonding partners V\ and v 2 for the corresponding inverse move and accepts 
the proposed move with probability: 



min(l,exp(-/5AE( J B / -»■ kA) + N b A/j,) 

B'^kA — 



PiGLCC 
R' 



v x v 2 

Thus the choice of the acceptance probabilities fullfills the detailed balance condition: 

P l A\r>P ro P pace _ p (r>l\-pprop pace 

r eq\^)^kA-^B ir kA^B' ~ r eq\D )^B'^kA r B'^kA 
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Figure 1: The measured line of pseudo phase coexistence separating the A-rich and A-poor 
phases, for which the concentration distribution has two peaks of equal weight. The position 
of the critical point as determined using the cumulant intersection method (see also figure 0) 
is indicated, as are the measured directions of the two relevant scaling fields h and r. 



Figure 2: The value of the fourth order cumulant ratio Gl = 1— < m 4 > /3 < m 2 > 2 with 
m = (p — 4> cx , expressed as a function of system size L and well depth e, along the line of pseudo 
phase coexistence. An intersection occurs for a value of Gl = 0.47 at e c = 0.02756(15), A/x c = 
0.003603(15) 



Figure 3: The normalised concentration distribution for the L = 40 and L = 64 system sizes 
obtained using the equal weight criterion with <fi* =< <fi >, at the assigned value of the critical 
well depth e = 0.02756(15). The data are expressed in terms of the scaling variable x = 
aJ^L^^^ — cpc), where the non-universal scale factor aj^ has been chosen so that distributions 
have unit variance. Also shown (solid) curve is the fixed point function p^ appropriate to the 
3D Ising universality class. Statistical errors do not exceed the symbol sizes. 



Figure 4: Extrapolation of (p cx (L), defined in the text, against Z^ 1 a >' u . The least squares fit 
yields an infinite volume estimate C = 0.03823(f9). 



Figure 5: The measured antisymmetric field mixing corrections 5pi{x) of the L = 40 and 
L = 64 critical concentration distributions, expressed in terms of the scaling variable x = 



a 



M 



\L^' v {(f) — 4> c ) and shown as the data points. The data has itself been corrected for a small 



off-coexistence correction as described in the text. Also shown (full curve) is the universal 
prediction following from equation |2.10| , utilising predetermined Ising forms 



Figure 6: The normalised distributions of the critical ordering operator Pl{M) for the L = 40 
and L = 64 system sizes, expressed as a function of the scaling variable y = aJ^L^^^Ai — -M c ). 
The full curve is the fixed point function p* M appropriate to the Ising universality class [p^j. 
The non-universal scale factor aj} has been chosen so that the distributions have unit variance. 
The data collapse corresponds to a choice of the field mixing parameter s = 0.06(1). Statistical 
errors do not exceed the symbol sizes. 



Figure 7: The normalised distributions of the critical energy density Pl{v) for the L = 40 and 
L = 64 system sizes. 



Figure 8: The normalised distribution of the energy-like operator Pl(£) (cf. equation |2.8b| ) 



expressed as a function of the scaling variable z = a £ 1 L^ 1 a ^ u (£—£ c ) for the L = 40 and L = 64 
system sizes. The full curve is the fixed point function pg appropriate to the Ising universality 



class |3ni. The non-universal scale factor a s , has been chosen so that the distributions have 



unit variance. The data collapse shown corresponds to a choice of the field mixing parameter 
r = —1.00(3). Statistical errors do not exceed the symbol sizes. 



Figure 9: The finite size scaling behavior of the variance of the critical energy, energy operator 
and ordering operator distributions, cf. figures [7] and [g. The straight lines superimposed 
on the data points have the forms 0.0076L 7 ^ (solid line) and 2.7L a ^ u (broken line), where 
7/1/= 1.970, a/u = 0.211. 



Figure 10: Illustration of the semi-grandcanonical Monte-Carlo moves for k = 3 and A?b = 3. 
Arrows indicate the bonds that are removed when one creates k = 3 A-chains or which have 
to be appropriately added in the inverse case. The monomer positions are left unaltered. 



